library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(reshape2)
library(ggplot2)
library(ggpmisc)# add linear regression equation to figures
## Loading required package: ggpp
## Registered S3 methods overwritten by 'ggpp':
## method from
## heightDetails.titleGrob ggplot2
## widthDetails.titleGrob ggplot2
##
## Attaching package: 'ggpp'
## The following object is masked from 'package:ggplot2':
##
## annotate
library(latex2exp)
library(ggpubr)# multiple figures in a page
##
## Attaching package: 'ggpubr'
## The following objects are masked from 'package:ggpp':
##
## as_npc, as_npcx, as_npcy
library(scales)
library(cowplot)# get_legend() & plot_grid() functions
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
##
## get_legend
N=10000
theta=7e-10 * 2e4 * 1e6
paras<-rbind(data.frame(Vs=1,sdMutation=0.05,speedCoef=c(0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.35,0.38,0.39,0.4,0.5,0.53,0.55,0.8,1.0,1.2,1.4,1.5,1.6)),
data.frame(Vs=1,sdMutation=0.1,speedCoef=c(0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.4,0.5,1.0,1.1,1.2,1.3,1.6,1.65,2.2,2.8,2.9,3.0,3.2,3.4,3.6)),
data.frame(Vs=1,sdMutation=0.5,speedCoef=c(0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.4,0.5,1.0,2.0,3.0,6.5,10.0,11.0,12.0,12.5,13.5,14.0,15.0,16.0,17.0,18.0,18.2,20.0,20.5,21.0)),
data.frame(Vs=10,sdMutation=0.05,speedCoef=c(0.002,0.005,0.01,0.02,0.05,0.07,0.08,0.09,0.1,0.11,0.2,0.3,0.4)),
data.frame(Vs=10,sdMutation=0.1,speedCoef=c(0.002,0.005,0.01,0.02,0.05,0.1,0.15,0.18,0.2,0.22,0.23,0.24,0.3,0.4,0.5,0.6,0.7,0.8)),
data.frame(Vs=10,sdMutation=0.5,speedCoef=c(0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.4,0.5,1.3,2.3,2.32,2.35,2.4,3.3,3.35,4.3,1.0,4.5,4.6,4.9,5.0,5.2,5.4,6.0)),
data.frame(Vs=100,sdMutation=0.05,speedCoef=c(0.002,0.005,0.01,0.02,0.025,0.03,0.031,0.035,0.04,0.05,0.1,0.2)),
data.frame(Vs=100,sdMutation=0.1,speedCoef=c(0.002,0.005,0.01,0.02,0.03,0.04,0.045,0.05,0.055,0.057,0.1,0.2,0.3)),
data.frame(Vs=100,sdMutation=0.5,speedCoef=c(0.002,0.005,0.01,0.02,0.05,0.1,0.2,0.3,0.35,0.38,0.39,0.4,0.5,0.52,0.55,0.7,0.9,1.1,1.3,1.4,1.5))
)
theme_set(theme_classic())
theme_update(text=element_text(size=13),
legend.text=element_text(size=10),legend.title=element_text(size=10,hjust=.5),
legend.key.size=unit(.3,'cm'),
legend.key.width=unit(.6,'cm'),
legend.spacing=unit(0,'cm'),
legend.background = element_rect(color="grey",fill=NA))
label_sdMutation<-function(string){TeX(paste("$\\omega : $",string))}
negativelog_trans <- function(base = exp(1)) {
trans <- function(x) log(-x, base)
inv <- function(x) -base^x
trans_new(paste0("reverselog_", format(base)), transform = trans, inverse=inv)
}
load data
load("data/processed/df2.RData")
Stationary state in the stabilizing-selection phase
# Lag is stationary
ggplot(SSE_df2,aes(x=generation/N-10,y=gap))+
geom_line()+
facet_wrap(sdMutation~Vs,scales="free_y",
labeller=labeller(Vs=label_both,sdMutation=as_labeller(label_sdMutation,label_parsed)))+
labs(x=TeX("$Generation\\ (\\times 10^4 )"),y="Lag")
Stationary state in the moving-optimum phase
# Lag is stationary
ggplot(MOE_df2,aes(x=generation/N-10,y=gap))+
geom_line(aes(group=speedCoef,col=speedCoef*sqrt(SSE_VG)/sdMutation))+
facet_wrap(sdMutation~Vs,scales="free_y",
labeller=labeller(Vs=label_both,sdMutation=as_labeller(label_sdMutation,label_parsed)))+
scale_color_viridis_c()+
labs(x=TeX("$Generation\\ (\\times 10^4 )"),y="Lag",col = "Speed")
Take the average of the 1000 generations
# take average
df<-group_by(MOE_df2,Vs,sdMutation,speedCoef)%>%
summarise(MOE_VG=mean(varPheno_mean),
MOE_VG_sd=sd(varPheno_mean),
MOE_gap=mean(gap),
MOE_gap_sd=sd(gap),
MOE_fitness=mean(meanFitness_mean))
## `summarise()` has grouped output by 'Vs', 'sdMutation'. You can override using
## the `.groups` argument.
df<-left_join(df,as.data.frame(parameters),by=c("Vs","sdMutation"))
df<-mutate(df,
speed=speedCoef * sqrt(SSE_VG),
sigmaOmega2=sdMutation^2/Vs,
group=paste0("Vs=",Vs,", $\\sigma_{m}=$",sdMutation))
df<-as.data.frame(df)
so=unique(df$sigmaOmega2)%>%sort
so=as.character(format(so,scientific = T))%>%unique
df$sigmaOmega2<-format(df$sigmaOmega2,scientific=T)
df$sigmaOmega2<-factor(df$sigmaOmega2,level=so,labels=so)
#critical speed at specific fitness cutoff
fit0.1<-filter(df,MOE_fitness>0.1)%>%group_by(Vs,sdMutation)%>%slice_max(speedCoef)%>%mutate(Cutoff="0.1")
fit0.01<-filter(df,MOE_fitness>0.01)%>%group_by(Vs,sdMutation)%>%slice_max(speedCoef)%>%mutate(Cutoff="0.01")
fit0.001<-filter(df,MOE_fitness>0.001)%>%group_by(Vs,sdMutation)%>%slice_max(speedCoef)%>%mutate(Cutoff="0.001")
fit0<-group_by(df,Vs,sdMutation)%>%slice_max(speedCoef)%>%mutate(Cutoff="0")
fit_df<-rbind(fit0.1,fit0.01,fit0.001)
# if simulate more, fit0.1$speedCoef<-c(0.38,1.15,12.2,0.072,0.185,2.35,0.022,0.042,0.37)
vcrit_p<-ggplot(fit_df,aes(x=sdMutation^2/Vs,y=speed/sdMutation/theta,col=Cutoff))+
geom_point()+
stat_poly_line()+
stat_poly_eq(use_label(c("eq", "R2")))+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10',breaks=c(0.01,0.1),labels=trans_format("log10", math_format(10^.x)))+
scale_color_manual(values=c("#F8766D","#00BA38","#619CFF"),breaks=c("0.001","0.01","0.1"),
labels=unname(TeX(c("$w_{th}>0.001","$w_{th}>0.01","$w_{th}>0.1"))))+
labs(x=TeX("$\\chi"),y=TeX("$Threshold\\ speed\\ (\\omega\\theta)"),col="Fitness threshold")+
theme(legend.position= c(0.7,0.2))+
annotation_logticks(sides = "bl")
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
vcrit_p_sigmaG<-ggplot(fit_df,aes(x=sdMutation^2/Vs,y=speedCoef))+
geom_point()+
stat_poly_line()+
stat_poly_eq(use_label(c("eq", "R2")))+
scale_x_continuous(trans='log10',label=comma)+
scale_y_continuous(trans='log10',label=comma)+
labs(x="Trait parameter",y=TeX("$Threshold\\ speed\\ (\\sigma_{G(S)})"))+
geom_hline(yintercept = c(0.04125416,0.79),linetype="dashed")+
theme(legend.position= c(0.7,0.2))
vcrit_p
vcrit_p_sigmaG
plot(fit0.1$speedCoef,fit0$speedCoef)
abline(0,1)
# pdf("output/threshold_speed_polygenic.pdf", width=4,height=3)
# vcrit_p
# dev.off()
ggplot(df,aes(x=speed,y=MOE_gap,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x="Shift speed",y="Lag")
Scaling by omega
lag_p<-ggplot(df,aes(x=speed/sdMutation/theta,y=MOE_gap/sdMutation,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y=TeX("$Mean\\ lag\\ (\\omega)"))
lag_p+theme(legend.position = "inside",legend.position.inside = c(0.2,0.6))
# pdf("output/lag_polygenic.pdf",width=4,height=3)
# lag_p+theme(legend.position = "inside",legend.position.inside = c(0.2,0.6))
# dev.off()
Scaling by sigma_g
lag_p2<-ggplot(df,aes(x=speed/sqrt(SSE_VG),y=MOE_gap,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y=TeX("$Lag"))
lag_p2
ggplot(df,aes(x=speed,y=MOE_fitness,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x="Shift speed",y='Fitness')
Scaling by omega
fit_p<-ggplot(df,aes(x=speed/sdMutation/theta,y=MOE_fitness,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y="Mean fitness")
fit_p+theme(legend.position = "inside",legend.position.inside = c(0.2,0.4))
# pdf("output/fitness_polygenic.pdf",width=4,height=3)
# fit_p+theme(legend.position = "inside",legend.position.inside = c(0.2,0.4))
# dev.off()
Scaling by sigma_g
fit_p2<-ggplot(df,aes(x=speed/sqrt(SSE_VG),y=MOE_fitness,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y="Mean fitness")
fit_p2
ggplot(df,aes(x=speed,y=MOE_VG,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x="Shift speed",y="Genetic variance")
Scaling by omega
vp_p<-ggplot(df,aes(x=speed/sdMutation/theta,y=MOE_VG/sdMutation^2,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y=TeX("$\\Genetic\\ variance\\ (\\omega^{2})"))
vp_p+theme(legend.position = "inside",legend.position.inside = c(0.8,0.4))
# pdf("output/Vg_polygenic.pdf",width=4,height=3)
# vp_p+theme(legend.position = "inside",legend.position.inside = c(0.8,0.3),
# legend.text=element_text(size=8),legend.title=element_text(size=8,hjust=.5),
# legend.key.size=unit(.2,'cm'),legend.key.width=unit(.4,'cm'))
# dev.off()
Scaling by sigma_g
vp_p2<-ggplot(df,aes(x=speed/sqrt(SSE_VG),y=MOE_VG,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y=TeX("$\\Genetic\\ variance"))
vp_p2
load data
load("data/processed/muta2.RData")
Take average of every generation at the last 1000 generations and the average of 100 replicates
moe_fixed_mean<-group_by(moe_fixed,Vs,sdMutation,speedCoef)%>%
summarise_at(vars(-replicate),mean,na.rm=T)
moe_fixed_sd<-group_by(moe_fixed,Vs,sdMutation,speedCoef)%>%
summarise_at(vars(-replicate),sd,na.rm=T)
colnames(moe_fixed_sd)[-c(1:3)]<-paste0(colnames(moe_fixed_sd)[-c(1:3)],"_sd")
moe_fixed2<-left_join(moe_fixed_mean,moe_fixed_sd,by=c("Vs","sdMutation","speedCoef"))
moe_fixed2<-left_join(moe_fixed2,parameters,by=c("Vs","sdMutation"))
moe_fixed2<-mutate(moe_fixed2,speed=speedCoef*sqrt(SSE_VG),
sigmaOmega2=sdMutation^2/Vs,group=paste0("Vs",Vs,"sd",sdMutation))
so=unique(moe_fixed2$sigmaOmega2)%>%sort
so=as.character(format(so,scientific = T))%>%unique
moe_fixed2$sigmaOmega2<-format(moe_fixed2$sigmaOmega2,scientific=T)
moe_fixed2$sigmaOmega2<-factor(moe_fixed2$sigmaOmega2,level=so,label=so)
df_moe_fixed2<-left_join(df,moe_fixed2[,1:(ncol(moe_fixed2)-4)],by=c("Vs","sdMutation","speedCoef"))
lag_effect_p<-ggplot(df_moe_fixed2,aes(x=MOE_gap/sdMutation,y=mean_effect/sdMutation,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10')+
labs(x=TeX("$Lag\\ (\\omega)"),y=TeX("$Mean\\ fixed\\ effect\\ size\\ (\\omega)"))
lag_effect_p2<-ggplot(df_moe_fixed2,aes(x=MOE_gap,y=mean_effect,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10')+
labs(x="Lag",y="Mean fixed effect size")
ggarrange(lag_effect_p,lag_effect_p2,ncol=1)
ggplot(moe_fixed2,aes(x=speed,y=n_fixed/1000,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x="Shift speed",y="Number of fixations")
# non-log x axis
ggplot(moe_fixed2,aes(x=speed,y=n_fixed/1000,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
labs(x="Shift speed",y="Number of fixations")
Scaling by omega
fn_p<-ggplot(moe_fixed2,aes(x=speed/sdMutation/theta,y=n_fixed/1000,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y="Number of fixations")
fn_p
# non-log x axis
fn_p_nonlog<-ggplot(moe_fixed2,aes(x=speed/sdMutation/theta,y=n_fixed/1000,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y="Number of fixations")
fn_p_nonlog
# pdf("output/fixation_number_polygenic.pdf",width=4,height=3)
# fn_p_nonlog+theme(legend.position = "inside",legend.position.inside = c(0.8,0.3))
# dev.off()
Scaling by sigma_g
fn_p2<-ggplot(moe_fixed2,aes(x=speed/sqrt(SSE_VG),y=n_fixed/1000,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y="Number of fixations")
fn_p2
# non-log x axis
ggplot(moe_fixed2,aes(x=speed/sqrt(SSE_VG),y=n_fixed/1000,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y="Number of fixations")
ggplot(moe_fixed2,aes(x=speed,y=mean_fixedtime,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x="Shift speed",y="Fixation time")
Scaling by omega
ft_p<-ggplot(moe_fixed2,aes(x=speed/sdMutation/theta,y=mean_fixedtime,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y="Fixation time")
ft_p
Scaling by sigma_g
ft_p2<-ggplot(moe_fixed2,aes(x=speed/sqrt(SSE_VG),y=mean_fixedtime,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y="Fixation time")
ft_p2
ggplot(moe_fixed2,aes(x=speed,y=mean_effect,group=group,col=sigmaOmega2))+
geom_hline(aes(yintercept=sdMutation * sqrt(2/pi)))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
annotate("text",x=0.8,y=0.05 * sqrt(2/pi)+0.005,label=TeX("$\\omega=0.05",output="character"),parse=T)+
annotate("text",x=0.8,y=0.1 * sqrt(2/pi)+0.02,label=TeX("$\\omega=0.1",output="character"),parse=T)+
annotate("text",x=0.8,y=0.5 * sqrt(2/pi)+0.05,label=TeX("$\\omega=0.5",output="character"),parse=T)+
labs(x="Shift speed",y=TeX("$Mean\\ \\alpha\\ of\\ fixations"))
Scaling by omega
mes_p<-ggplot(moe_fixed2,aes(x=speed/sdMutation/theta,y=mean_effect/sdMutation,group=group,col=sigmaOmega2))+
# geom_hline(yintercept=sqrt(2/pi))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
#scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
#scale_y_continuous(trans='log10')+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y=TeX("$Mean\\ \\alpha\\ of\\ fixations\\ (\\omega)"))
mes_p
# pdf("output/fixation_effect_polygenic.pdf",width=4,height=3)
# mes_p+theme(legend.position = "inside",legend.position.inside = c(0.8,0.3))
# dev.off()
Scaling by sigma_g
mes_p2<-ggplot(moe_fixed2,aes(x=speed/sqrt(SSE_VG),y=mean_effect,group=group,col=sigmaOmega2))+
# geom_hline(yintercept=sqrt(2/pi))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y=TeX("$Mean\\ \\alpha\\ of\\ fixations"))
mes_p2
ggplot(moe_fixed2,aes(x=speed,y=sd_effect,group=group,col=sigmaOmega2))+
geom_hline(aes(yintercept=sdMutation * sqrt(1-2/pi)))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
annotate("text",x=0.8,y=0.05 * sqrt(1-2/pi)+0.005,label=TeX("$\\omega=0.05",output="character"),parse=T)+
annotate("text",x=0.8,y=0.1 * sqrt(1-2/pi)+0.01,label=TeX("$\\omega=0.1",output="character"),parse=T)+
annotate("text",x=0.8,y=0.5 * sqrt(1-2/pi)+0.02,label=TeX("$\\omega=0.5",output="character"),parse=T)+
labs(x="Shift speed",y=TeX("$Sd(\\alpha)\\ of\\ fixations"))
Scaling by omega
sde_p<-ggplot(moe_fixed2,aes(x=speed/sdMutation/theta,y=sd_effect/sdMutation,group=group,col=sigmaOmega2))+
# geom_hline(yintercept=sqrt(1-2/pi))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y=TeX("$Sd(\\alpha)\\ of\\ fixations\\ (\\omega)"))
sde_p
Scaling by sigma_g
sde_p2<-ggplot(moe_fixed2,aes(x=speed/sqrt(SSE_VG),y=sd_effect,group=group,col=sigmaOmega2))+
# geom_hline(yintercept=sqrt(1-2/pi))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y=TeX("$Sd(\\alpha)\\ of\\ fixations"))
sde_p2
ggplot(moe_fixed2,aes(x=speed,y=mean_posEffect,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
labs(x="Shift speed",y="Mean positive fixed effect")
Scaling by omega
ggplot(moe_fixed2,aes(x=speed/sdMutation/theta,y=mean_posEffect/sdMutation,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y=TeX("$Mean\\ positive\\ fixed\\ effect\\ (\\omega)"))
ggplot(moe_fixed2,aes(x=speed,y=mean_negEffect,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans=negativelog_trans(10),breaks=c(0,-0.1,-0.2,-0.4,-0.6))+
labs(x="Shift speed",y="Mean negative fixed effect")
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
Scaling by omega
ggplot(moe_fixed2,aes(x=speed/sdMutation/theta,y=mean_negEffect/sdMutation,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans=negativelog_trans(10),breaks=c(0,-0.25,-0.75,-1))+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y=TeX("$Mean\\ negative\\ fixed\\ effect\\ (\\omega)"))
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
Take average of every generation at the last 1000 generations and the average of 100 replicates
moe_muta_mean<-subset(shift_muta2,generation>19.9*N)%>%
group_by(Vs,sdMutation,speedCoef,replicate)%>%
summarise_at(vars(-generation),mean,na.rm=T)%>%
group_by(Vs,sdMutation,speedCoef)%>%summarise_at(vars(-replicate),mean,na.rm=T)
moe_muta_sd<-subset(shift_muta2,generation>19.9*N)%>%
group_by(Vs,sdMutation,speedCoef,replicate)%>%
summarise_at(vars(-generation),mean,na.rm=T)%>%
group_by(Vs,sdMutation,speedCoef)%>%summarise_at(vars(-replicate),sd,na.rm=T)
colnames(moe_muta_sd)[-c(1:3)]<-paste0(colnames(moe_muta_sd)[-c(1:3)],'_sd')
moe_muta2<-left_join(moe_muta_mean,moe_muta_sd,by=c('Vs','sdMutation','speedCoef'))
moe_muta2<-left_join(moe_muta2,parameters,by=c('Vs','sdMutation'))
moe_muta2<-mutate(moe_muta2,speed=speedCoef * sqrt(SSE_VG),
sigmaOmega2=sdMutation^2/Vs,
group=paste0("Vs=",Vs,", $\\sigma_{m}=$",sdMutation))
so=unique(moe_muta2$sigmaOmega2)%>%sort
so=as.character(format(so,scientific = T))%>%unique
moe_muta2$sigmaOmega2<-format(moe_muta2$sigmaOmega2,scientific=T)
moe_muta2$sigmaOmega2<-factor(moe_muta2$sigmaOmega2,level=so,label=so)
ggplot(moe_muta2,aes(x=speed,y=n_segregations,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x="Shift speed",y="Number of segregating sites")
Scaling by omega
sn_p<-ggplot(moe_muta2,aes(x=speed/sdMutation/theta,y=n_segregations,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y="Number of segregating sites")
sn_p
# pdf("output/segregating_number_polygenic.pdf",width=4,height=3)
# sn_p+theme(legend.position = "inside",legend.position.inside = c(0.8,0.8),
# legend.text=element_text(size=8),legend.title=element_text(size=8,hjust=.5),
# legend.key.size=unit(.2,'cm'),legend.key.width=unit(.4,'cm'))
# dev.off()
Scaling by sigma_g
sn_p2<-ggplot(moe_muta2,aes(x=speed/sqrt(SSE_VG),y=n_segregations,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y="Number of segregating sites")
sn_p2
ggplot(moe_muta2,aes(x=speed,y=n_posSegregations/n_segregations,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(labels=percent)+
labs(x="Shift speed",y="Proportion of positive segregating mutations")
Scaling by omega
ggplot(moe_muta2,aes(x=speed/sdMutation/theta,y=n_posSegregations/n_segregations,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(labels=percent)+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y="Proportion of positive segregating mutations")
Scaling by sigma_g
ggplot(moe_muta2,aes(x=speed/sqrt(SSE_VG),y=n_posSegregations/n_segregations,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(labels=percent)+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y="Proportion of positive segregating mutations")
ggplot(moe_muta2,aes(x=speed,y=mean_effect,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
labs(x="Shift speed",y=TeX("$Mean\\ \\alpha\\ of\\ segregating\\ sites"))
Scaling by omega
mse_p<-ggplot(moe_muta2,aes(x=speed/sdMutation/theta,y=mean_effect/sdMutation,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y=TeX("$Mean\\ \\alpha\\ of\\ segregating\\ sites\\ (\\omega)"))
mse_p
# pdf("output/segregating_effect_polygenic.pdf",width=4,height=3)
# mse_p+theme(legend.position = "inside",legend.position.inside = c(0.8,0.35),
# legend.text=element_text(size=8),legend.title=element_text(size=8,hjust=.5),
# legend.key.size=unit(.2,'cm'),legend.key.width=unit(.4,'cm'))
# dev.off()
Scaling by sigma_g
mse_p2<-ggplot(moe_muta2,aes(x=speed/sqrt(SSE_VG),y=mean_effect,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y=TeX("$Mean\\ \\alpha\\ of\\ segregating\\ sites"))
mse_p2
ggplot(moe_muta2,aes(x=speed,y=mean_posEffect,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
labs(x="Shift speed",y="Mean positive segregating effect")
Scaled by omega
ggplot(moe_muta2,aes(x=speed/sdMutation/theta,y=mean_posEffect/sdMutation,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y=TeX("$Mean\\ positive\\ segregating\\ effect\\ (\\omega)"))
Scaled by sigma_g
ggplot(moe_muta2,aes(x=speed/sqrt(SSE_VG),y=mean_posEffect,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans='log10')+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y=TeX("$Mean\\ positive\\ segregating\\ effect"))
ggplot(moe_muta2,aes(x=speed,y=mean_negEffect,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans=negativelog_trans(10))+
labs(x="Shift speed",y="Mean negative segregating effect")
Scaled by omega
ggplot(moe_muta2,aes(x=speed/sdMutation/theta,y=mean_negEffect/sdMutation,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans=negativelog_trans(10))+
labs(x=TeX("$Shift\\ speed\\ (\\omega\\theta / generation)"),y=TeX("$Mean\\ negative\\ segregating\\ effect\\ (\\omega)"))
Scaled by sigma_g
ggplot(moe_muta2,aes(x=speed/sqrt(SSE_VG),y=mean_negEffect,group=group,col=sigmaOmega2))+
geom_point()+
geom_line()+
scale_color_viridis_d(name=TeX("$\\chi"),begin=0,end=0.9)+
scale_x_continuous(trans='log10',labels=trans_format("log10", math_format(10^.x)))+
scale_y_continuous(trans=negativelog_trans(10))+
labs(x=TeX("$Shift\\ speed\\ (\\sigma_{G(S)} / generation)"),y="Mean negative segregating effect")